function [W1, W2] = BackpropBatchCE(W1, W2, X, D)
  alpha = 0.9;
  
  N = 4;  
  dW1sum = zeros(size(W1));
  dW2sum = zeros(size(W2));
  for k = 1:N
    x = X(k, :)';        % x = a column vector
    d = D(k);
    
    v1 = W1*x;
    y1 = Sigmoid(v1);    
    v  = W2*y1;
    y  = Sigmoid(v);
    
    e     = d - y;
    delta = e; % CE

    e1     = W2'*delta;
    delta1 = y1.*(1-y1).*e1; 
    
    dW1sum = dW1sum + alpha*delta1*x';
    
    dW2sum = dW2sum + alpha*delta*y1';    
  end
  dWavg1 = dW1sum / N;
  dWavg2 = dW2sum / N;
  W1 = W1 + dWavg1;
  W2 = W2 + dWavg2;
end